function components = obtain_wav_frequency_components(wav_fullname)
% Algorithm that calculates the main frequency components from 
% a wav filename and returns an array with the frequency and the 
% amplitude.
%
% components = [ frequency 1, amplitude 1; frequency 2, amplitude 2; ...]
%


%% Read a wav file
%Fs: Sampling frequency, nbits:number of bits per sample
[y, Fs, nbits] = wavread(wav_fullname);

T = 1/Fs;            % Sample time
L = length(y);      % Length of signal
%t = (0:L-1)*T;      % Time vector

%% Select interval

interval_duration = 0.5; %seconds
number_of_samples = round(interval_duration/T);
start_sample = round(L/2);
yy = y(start_sample: start_sample+number_of_samples-1);
LL = length(yy);

%% FFT 
%NFFT = 512; 
NFFT = 2^nextpow2(LL);
Y = fft(yy,NFFT)/LL;
f = Fs/2*linspace(0,1,NFFT/2);

%plot(f,2*abs(Y(1:NFFT/2)));

%% Obtain peaks

npoints = NFFT/2;
[val,pos,YY] = obtain_peaks(Y,npoints);
peak_freq = round(f(pos));
%disp(['First peak:' num2str(val) ' ' num2str(peak_freq) ' Hz']);
components = [peak_freq val];

for i=2:5
[val,pos,YY] = obtain_peaks(YY,npoints);
peak_freq = round(f(pos));
components = [components; peak_freq val];
end

% for i=1:size(components,1)
%     disp(['Peak:' num2str(components(i,2)) ' ' num2str(components(i,1)) ' Hz']);
% end

end 


function [val,pos,YY] = obtain_peaks(Y,npoints)

YY = Y;
[val,pos]=max(abs(YY(1:npoints)));
%disp(['Peak:' num2str(val) ' ' num2str(peak_freq) ' Hz']);
peak_amplitude = 10;

p1 = pos - peak_amplitude;
p2 = pos + peak_amplitude;
if p1 < 0
    p1 = 1;
end
if p2 > npoints
    p2 = npoints;
end
for i=p1:p2
    YY(i) = 0;
end

end
